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1. INTRODUCTION 


The Lockheed Missiles and Space Company's SL4 Code is a computer implemented 
computational scheme which provides 'theoretical predictions of the radiative 
and convective energy transfer to the surface of a blunt body during hyper- 
velocity entry into planetary atmospheres. Although prediction of the convective 
heating rate is a direct consequence of the solution method, the SL4 Code was 
developed primarily to analyze those entry conditions where radiative processes 
dominate the total energy transport. The surface heating rate prediction 
results from a solution of the coupled equations describing the flow of a 
viscous, radiating gas past a blunt body which may or may not be undergoing 
massive ablation. The detailed derivation and analysis of those equations can 
be found in Ref. 1. 

The solution provided by the SL4 Code is a solution to the inverse problem. 
In order to use this code, the shock wave must be specified completely. The 
body which supports this shock will be a part of the solution.. In the SL4 
Code, a specific shock equation was built-in. This equation contains two 
free parameters. Detailed discussion on this equation will be given in 
Section 3 . It is believed that these two families of shock waves will be 
sufficient to describe many cases of practical interest. 

The SL4 Code also requires an a priori specification of the streamwise 
variation of the surface blowing rate. A particular equation for the blowing 
rate distribution, which contains one free parameter, was built into the code. 
This equation will be discussed in Section 3- By performing an iterative 
process, the blowing rate distribution can be made to reflect the surface 
radiative flux. 

Within the context of its application by a knowledgeable user, the 
code is operational on the UNIVAC 1108 (Exec. 8 ) system in use at IMSC and the 
CDC 76 OO system in use at the University of California, Berkeley, California. 
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2. COMPUTATION STRUCTURE 


A chart is given in Fig. 1 which shows the logic sequence followed in 
obtaining an iterative solution to the conservation equations. Details on the 
input/output quantities will be given in the next section. The first step 
in the SL4 Code computation is the determination of the shock layer parameters 
p's. For the definitions of P's, the reader should consult Ref. 1. The 
P's are determined by the two separate subroutines. One subroutine supplies 
the shock shape and its derivatives, while the other subroutine describes the 
injection rate and computes the total injected mass and its derivatives. The 
next step in the computation is the determination of the properties behind the 
shock. From the oblique shock relations, the pressure and the tangential 
velocity are determined. After neglecting the normal velocity and the radiative 
energy loss right behind the shock, the enthalpy can be determined. A simple 
iteration using the FEMP subroutine to calculate the actual post shock density 
leads rapidly to a determination of all shock properties. 

The sequential solution of the conservation equations is now initiated. 

The input estimates of the tangential velocity, blown gas mass fraction, total 
enthalpy, and their streamwise derivatives provide the data required to start 
the iteration process. The first step in the iterative cycle is to calculate 
collisional transport properties. This is accomplished using the FEMP package 
set of subroutines. Along with the enthalpy and elemental mass fractions (as 
deduced from the local value of the blown gas mass fration), the FEMP routine 
requires a value for the pressure. This value is given by an analytic expression 
as presented in Ref. 1. From FEMP, the variation in temperature, density, 
viscosity and binary diffusion coefficient across the shock layer is obtained. 
These properties provide the data required to perform a solution to the momentum 
equations using the technique discussed in Ref. 1. Convergence of the velocity 
profile is obtained via a local iterative loop. For a fixed set of density and 
diffusion coefficient data together with the just -calculated velocity distribution, 
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Fig. 1. Computational Logic Flow 
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the diffusion equation can be solved in a similar manner. 


At this point our procedure elects to update the species number density- 
distribution to reflect the change in the blown gas mass fraction profile. It 
is essential to recognize that in flows dominated by radiative transport 
effects, temperature is a crucial variable. Indeed, the temperature which 
is calculated in FEMP from the input enthalpy-mass fraction distribution must 
be a "reasonable" distribution. Having established a reasonable temperature 
distribution from the initial guess it is imperative to the stability of the 
solution to prevent the initial temperature profile from being altered to a 
large degree in a single step. Such a large temperature change often will 
occur if we hold the enthalpy at its initial value and calculate a new temperature, 
density and species number density profile using the updated elemental mass 
fraction distribution. The reason which underlies the temperature change problem 
is the following: For flows dominated by radiative transport and especially at 

high Reynolds numbers, the shear layer is nearly isothermal. However, a rapid 
change in elemental mass fraction occurs across the shear layer. Since the 
effective heat capacity of the blown gas may be quite different than that of the 
atmospheric gas, then the rapid change in the elemental mass fractions is 
accompanied by rapid change in the enthalpy despite the nearly isothermal 
character of the shear layer. Hence, if in the updating of the elemental mass 
fractions, the injected gas moves to a point where the enthalpy is significantly 
different, the resultant temperature is physically unrealistic and causes 
catastrophic failure of the energy equation. 

To eliminate this high sensitivity between the diffusion and energy 
equations, our approach fixes temperature at the current value obtained on the 
first pass through FEMP. Thus in the second pass through FEMP, the temperature 
is the independent thermodynamic variable and updated values of enthalpy and 
species number density are determined. Of course, for a nonblowing there is 
no requirement for a second pass through FEMP. Thus, for a nonblowing problem, 
these two stages are bypassed as indicated in Fig. 1. 

Having updated the species number densities and using the current 
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temperature profile, the next step is to determine the radiative flux diver- 
gence across the shock layer. In deriving the equations describing the 
radiative transport across the shock layer, it was convenient to use the 
physical shock layer normal coordinate. Accordingly, the first task in the 
radiative transport subroutine is to determine from the m coordinate -values 
and the density distribution the physical thickness of the shock layer. The 
radiative flux divergence is then calculated at every uj point. These flux 
divergence values are used along with the current velocity profile and Prandtl 
number variation to provide a new estimate for the enthalpy profile. 

A test for convergence of the total enthalpy profile is made at this point 
in the computational cycle. Further comments on the level of convergence to 
be expected will be given in Section 4. If satisfactory convergence has not 
been achieved, the entire iterative procedure is reinitiated. A number of 
investigators besides the author have observed the instability which arises 
if the calculated enthalpy profile is used unaltered as the enthalpy profile 
for the subsequent iterative cycle. In the absence of a meaningful stability 
analysis (which is likely to be impossible to achieve in view of the strong 
nonlinearity of the radiative transport terms), the following expedient has 
been employed. Let H 1 be the current guess on the enthalpy profile and H i+1 
the resulting calculated value. Then the next guess on the enthalpy profile 
which will be used to reinitiate the iterative loop is 

H = (l-m)H 1 + mH 1+1 
next 

where m is a fraction less than 1.0. A practical upper limit on m appears 
to be 0. 5- As m -» 0 stability will undoubtedly be achieved but at the cost 
of numerous iterations for significant progress toward convergence. Our ex- 
perience has shown that m ~ 0.2 to 0.3 are reasonable values for blowing problems, 
while m may approach 0.5 for nonblowing. 

If convergence is achieved or, as is more likely, a maximum number of 
iterative cycles have passed, then a final pass through the radiative transport 
subroutine is made. In this final pass the radiative flux at the wall is 
calculated using the unaltered enthalpy profile, H = H i+ \ This final wall 
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radiative flux value then can be compared against the value obtained from 
the enthalpy profile which existed upon the last entry to the energy equation. 

The comparison of these two flux values provides the best measure of convergence. 
Preparation of data for output, which is discussed in Section 3> completes 
the computation. 
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3. INPUT AND OUTPUT 


3 • 1 Input Data 

Table 1 lists the input data required to initiate a SL4 Code run. Also 
given in Table 1 is the data format, FORTRAN symbol and a synoptic description 
of each input quantity. The data is segregated into card groups. Depending 
on the amount of data required, more than one punched card may be needed in a 
given group. Figures 2 and 3 are copies from a computer listing of a set of 
input data. In the following discussion we shall proceed sequentially through 
Table 1. The reader will find it instructive to examine the typical values 
shown in Figs. 2 and 3 as we move through Table 1. The centimeter-gram-second 
system of units, along with °K for temperature, has been used for those input 
quantities given in dimensional form. 

The data specified in the first card group is self-explanatory. The data 
specified in the second card group deals primarily with description of the 
free stream properties. Vehicle velocity and ambient density are self-explana- 
tory. The shock bluntness p and the asymptotic shock angle 9 are the two 

s s 

free parameters in describing the shock shape. In the SlA Code, this shock 
shape is given by the following equation. 

2 2 
r = 2z + r z tan 9 - B z 
s s s s 

Where z is the horizontal distance from shock stagnation point, and r^ is 

the normal distance from the axis of symmetry. For 9 =0, this equation 

s 

describes: a parabola if B = 0; a sphere if 3 = 1; a hyperbola if P <0; 

s s s 

and an ellipsoid if 3 >1. For a nonzero value of 9 , if 3 =0, .the shock shape 

s s s 

equation describes a parabola near the nose and then asymptotically the shock 
shape approaches a cone with half angle equal to 9 g at large z . 

The wall is characterized by either a prescribed enthalpy or (normally) 
temperature according to an option given in card group 6. If a prescribed 
wall enthalpy is input, that enthalpy must be normalized to the total free 
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2 

stream enthalpy + u <J2. If a prescribed wall temperature is used, the 
value is input in dimensional (°K) form. 


The blowing rate is given in terms of a normalized value m = p v /p u 

V W oo oo 

which is given by the following equation as a function of r 

s 


m 


m 

o 


1+Ar 2 

s 


where m Q is the blowing rate at stagnation point and the parameter A 
characterizes the rate change of the local blowing rate around the body. 

The convergence criterion t is applied by requiring 
H 1+1 - H 1 , r 
(H 1+1 + ^)/2 

be satisfied at all shock layer points in order for convergence to be achieved. 
For nonblowing problems, convergence to a level of T = .01 (i.e., 1$ convergence) 
is an achievable goal. However, for radiative dominated problems, a 1# 
convergence tolerance is probably too severe. Usually a few points in the 
shock layer hold -up convergence. Hence, when a 1$ convergence level is achieved 
the radiative flux at the wall is changing considerably less than 1% between 
iterations. A precision of 1$ is inconsistent with the overall accuracy of 
the optical properties, transport properties and the assumption of chemical 
equilibrium which set the true error bound on the heating predictions. For 
blowing problems, convergence to a 5$ level is rarely achieved. For blowing 
problems, convergence is much more reasonably measured by comparing the wall 
flux based on the final guessed and calculated enthalpy profile. 


The weighting factor employed in stabilizing the successive iterations 
on the energy equation was discussed in the previous section. 

The FEMP routine considers two different temperature ranges. These 
two ranges are denoted as the low temperature regime and the high temperature 
regime. The low temperature regime does not contain ionized species whereas 
the high temperature regime does not contain polyatomic species. A tabulation 
of the species considered in each temperature regime is given in Table I of 
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Table 1 

SL4 Code Input Requirement 


Card Group 

Format 

Item 

1 

6a6 

RUNID 

2 

7E10.4 

UIN 


RIN 

RS 


BS 

THS 

THW 

XMZ 

3 4E10.4 AAA 

CORR 

CONV 

TLH 


4 5E10.4 ALPV(l) 

ALPV(2) 
ALPV(3) 
ALPV(4) 
ALPV( 5) 

5 5E10.4 ALPA(l) 

ALPA(2) 

ALPA(3) 

ALPA(4) 

alpa( 5 ) 

6 715 NRS 

ITGS 

IDG 

ICWT 

MOL 


Description 

Alphanumeric run identification 

Vehicle velocity (cm/sec) 

Ambient density (gm/cm3) 

Shock radius of curvature at 
stagnation point 
Shock bluntness 
Asymptotic shock angle (deg) 

Wall enthalpy or wall temperature ( °K) 
(see ICWT, card 6) 

Blowing rate at stagnation point 

Blowing rate parameter 
Weighting factor in enthalpy 
profile guess 

Convergence criteria enthalpy 
profile 

Low to high species temperature 
switch (°K) 

Elemental mass fraction of H at wall 
Elemental mass fraction of He at wall 
Elemental mass fraction of C at wall 
Elemental mass fraction of N at wall 
Elemental mass fraction of 0 at wall 

Elemental mass fraction of H at shock 
Elemental mass fraction of He at shock 
Elemental mass fraction of C at shock 
Elemental mass fraction of N at shock 
Elemental mass fraction of 0 at shock 

Number of z points input (max. 10 ) 
Maximum number of iterations 
= 0 no diagnostic output; =1 diag- 
nostic output 
= 0 h input; =1 T input 
= 0 without molecules; =1 with molecules 
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Table 1 (Continued) 


"d Group 

Format 

Item 

Description 

6 

715 

LINES 

IRAD 

= 0 without lines ;=1 with lines 
= 0 no radiation; =1 with radiation 

7 

8E10.4 

ZSS(l),I=l,NRS 

z points 

8 

15 

NES 

number of grid points input 

9 

8eio.4 

a>(K),K=l,NES 

uu grid points 

10 

Seio.4 

sf(k),k=i,nes 

Velocity profile~guess 

11 

8E10.4 

c(k),k=i,nes 

Blown gas mass fraction' profile 
guess (if XMZ=0, input C(K)=0) 

12 

8eio. 4 

h(k),k=i,nes 

Enthalpy profile guess 

13* 

8E10.4 

eta(k),k=i,nes 

Stream function guess 

14* 

8E10.4 

cf(k),k=i,nes 

F 

15* 

8eio.4 

cx(k),k=i,nes 

C 

V 

1 6 * 

8E10.4 

hx(k),k=i,nes 

H 


* 

Input only if 


ZSS(l) > 10 
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Fig. 2. Input Data Vehicle and Atmospheric Parameters 
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Fig. 3. Input Data: Initial Shock Layer Profiles 



Ref. 2. Changes to that tabulation since it was published include: (l) addition 

of C^, He, HgO and CO^ to the low temperature regime; (2) addition of C^, He and 
He + to the high temperature regime. The SL4 code user must specify a tempera- 
ture at which the FEMP calculation is switched from the low to high regimes. 
Selection of this low to high temperature break-point must take into considera- 
tion that it should be sufficiently high that polyatomic species would not be 
present in significant amounts (< 1$) at or above that temperature nor ionized 
species present in significant amounts at or below that temperature. The 
temperature break-point will depend on the stagnation-point pressure and 
typically will vary from 5000°K at low pressures (~ 0.1 atm) to 8000°K at high 
pressures (~ 100 atm) . The solution is not critically dependent on the value 
of the break-point temperature selected. With a little experience in a given 
problem, an appropriate value for the temperature switch will be readily 
apparent to the user. An additional requirement is that the low to high 
temperature switch be greater than the wall temperature but less than the 
shock temperature. 

The fourth and fifth card groups permit specification of an arbitrary 

ablation product (card group 4) and atmospheric (card group 5) gas composition. 

A fundamental restriction on the SL4 code is that the composition of ablator 

and atmosphere be comprised of H, He, C, N and 0 elements. The species included 

in the FEMP routine were selected on the basis of being the most significant 

2 3 

for a mixture comprised of these elements. For both ablator or atmosphere 
the composition is given simply on an elemental mass fraction basis. The 
actual species composition of the gases entering the shock layer is not specified 
as this will be handled by FEMP. 

The sixth card group provides control options. The number of shock layer 
points is at the discretion of the user. Storage limitations require that the 
number of shock layer points be kept equal to or less than 40. The maximum 
number of iterations provides a termination point for the calculation when 
"convergence" of the enthalpy profile is not achieved. On the UNIVAC 1108 
Exec 8 system in operation at IMSC, the computation time for a full iterative 
cycle (i.e., solution of momentum, mass and energy equations including FEMP 
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and TRANS subroutines) for a typical problem involving 26 shock layer points 
is approximately 2 minutes. Computer time, then, is normally the overriding 
consideration in selection of the maximum number of iterations. 

The remaining quantities in card group 6 are self-explanatory options. 

With regard to the diagnostic output, there will not be a discussion of the 
diagnostic data in Section 3-2. Without extensive study of the code details, 
the diagnostic data would have no utility to the user. 

The last sequence of card groups specifies the computational grid 
(in the uj- coordinate ) and initial estimates of the normalized distribution 
across the shock layer of the stream function 7), tangential velocity f and 
its streamwise derivative F , the blowing gas mass fraction and its 
streamwise derivative , the total enthalpy H and its streamwise derivative 
H. Near the stagnation point (z < 10 ) no estimates will be necessary for 
T], F, and H. In the stagnation region F, C y and H are zero. These estimates, 

of course, are crucial and a discussion of input profiles will be given in 
Section 4. 

In addition to the normal input listed in Table 1, the user may also 
specify the spectral emissivity of the wall to account for surface emission 
and reflection. It was felt that it would be an inconvenience to require this 
spectral emissivity data to be input for every run. In a typical application 
of the SL4 code, we expect the surface properties to be fixed for a number of 
cases in which the atmospheric and vehicle parameters are varied. Accordingly, 
the spectral emissivity data are specified in an internally programmed data 
statement at the beginning of TRANS, the radiative transport subroutine. The 
data is specified in terms of the average emissivity within each of the 26 
spectral groups employed in the continuum transport evaluation. That is to say, 
the SL4 code currently represents the spectral variation of surface emissivity 
by average values at 2 6 locations across the spectrum. The frequency values of 
the continuum groups are evident from the output discussed in the next subsection. 
For the line transport, the SL4 code lifts the surface emissivity values appro- 
priate to the frequency interval of the various line groups in exactly the 
same way it determines the appropriate underlying continuum optical depth for 
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each line group. Currently, the SL4 code is programmed with unit surface 
emissivity for each continuum spectral interval. 

3 • 2 Output Data 

Figures 4 through 13 display copies of the computer listing of the 
output data for the problem defined by the set of input data shown in Figs. 2 
and 3* In Fig. 4 the thermodynamic properties behind the normal shock are shown. 
In addition, the shock layer thickness, the normalized tangential velocity, 
the shock, body and the interface coordinate are output. Also included on 
Fig. 4 is a status report on the last energy equation iteration in terms of a 
comparison between the guessed and calculated enthalpy profile. 

Figure 5 summarizes the radiative flux picture. Both the continuum and 
2 

line fluxes (in W/cm ) at the wall Q + (0) and shock Q - (l) are given as a 
function of the continuum and line group and, for convenience, the average 
spectral values (in W/cm eV) in each group were output. The final quantity 
in Fig. 5 is the integral of the radiative flux divergence over the shock layer. 
This integral is to be compared with the sum of the continuum and line total 
fluxes (i.e., spectral summation) emergent from shock layer at the wall and 
shock. Agreement of the flux divergence integral with the total summation of 
all flux components demonstrates adequate spatial resolution of the flux 
divergence calculation across the shock layer. 

The radiative flux is calculated at the same time the flux divergence 

is calculated. This allows us to present the radiative flux distribution 

across the shock layer in terms of the flux values at each u)-point. Figure 6 

shows a representative set of radiative flux data. The <u coordinate and 

corresponding physical coordinate (normalized to total shock layer thickness) 

are indicated. This is followed by a listing of the continuum and line fluxes, 

2 

in W/cm , and the spectral summation. The spectral continuum optical depth 
evaluated at the midpoint of each continuum or line group is presented along 
with the flux values. 
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At the termination of the radiative flux output at the last' shock layer point, 

TV = 1, the wall convective . and radiative heating values are presented. This is shown 
in Fig. 7- Of course, the radiative flux value was already output in terms 
of the continuum and line components in Fig. 5 . These two radiative flux 
components are simply summed and repeated for convenience in the output 
indicated by Fig. 7- In addition to the absolute heating values, a convective 
and radiative heat transfer coefficient, defined as the heating rates normalized 
by 1/2 P<x> u ® > : '' s P resen f e< I as shown in Fig. 7- 

A summary of the final values of the gasdynamic quantities--stream 
function, tangential velocity, enthalpy, ablation product mass fraction, 
flux divergence (w/cnr) is output next. A representative set is shown as 
Fig. 8 . A quick comparison between the initial and final profiles is thus 
available by examining the data of Figs. 3 and 8 . It will be observed that 
the enthalpy tabulated in Fig. 8 is the weighted average of the final guessed 
and calculated values shown in Fig. 4. Thus this enthalpy profile represents 
the input to the next iteration. In principal, the output shown in Fig. 8 
can be used directly to initiate a new set of iterations. 

A summary of the thermodynamic and collisional transport properties 
consistent with the gasdynamic data given in Fig. 8 is the next data output. 

A representative set is shown as Fig. 9- 

The next series of output presents the molecular and atomic species number 
density distributions across the shock layer. A representative set for the 
molecular species is shown as Fig. 10 and likewise for the atomic species as 
Fig. 11. 


The final computation which the SL4 code performs is an evaluation of the 
radiative flux at the wall using the unaltered calculated enthalpy profile 
(as opposed to the profile weighted between guessed and calculated values) 
from the last iteration. The representative values of output which describes 
this final calculation is shown as Fig. 12. In Fig. 12, data from the FEMP 
routine is displayed with the temperature variable being of maximum interest. 
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The final output gives the continuum and line flux which may be compared with 
the data shown on Fig. 5 to assess the level of convergence of the overall 
radiative flux determination. 
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Fig. 4. Output Data: Shock Conditions and Status of Enthalpy Iteration 
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Figure 5 

Output Data: Radiative Flux at Wall and Shock 
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Output Data: Typical Radiative Flux Data at Selected Shock Layer Points 
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Output Data: Shock Layer Gasdynamic Variables 
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Output Data: Shock Layer Thermodynamic and Transport Variables 
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Figure 10 

Output Data: Shock Layer Polyatomic and Diatomic Species 
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Figure 11 

Output Data: Shock Layer Atomic Species 
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Output Data: Final Temperature /Enthalpy Profile and 

Resultant Surface Flux 



it-. PRACTICAL ADVICE TO THE USER 


Set forth in this section are some practical guides to using the SL4 Code 
as deduced by the author as a result of many attempts at obtaining heating 
predictions . 

The fundamental problem is prescribing a set of initial profiles 
(velocity, blown gas mass fraction, enthalpy) in a nonphysical coordinate 
which is sufficiently close to reality to get underway with a series of 
iterations. In this regard, one should start the code at stagnation point. 

At stagnation point, only f, and H profiles need to be guessed. In 
Ref. k a detailed discussion was made on how to guess the profiles effectively 
at the stagnation point. Same rules should apply to the use of the SLb 

g 

Code. The term stagnation point is defined as z = 10 in this code, since, 
when z is exactly zero, the shock layer parameters. P's, are all undetermined. 

Once the stagnation point solution is obtained, one can move a finite 
distance away from the stagnation point, and use the stagnation point solution 
as the initial guessed profiles. In this manner the user can build-up 
solutions around the body. , The distance one moves from point to point in 
obtaining the local solutions will be dependent on the shock shape and the 
blowing variation. In general one should move only a small distance near the 
stagnation region. The size of the distance between solutions, Az, should be 
inversely proportional to the shock bluntness. That is, for more blunted 
shocks, Az should be smaller than the value for a less blunted shock. Also, 

Az is inversely proportional to the rate change of mass injection rate. 

The faster the injection rate changes, a smaller Az must be chosen. From the 
sample calculations of Ref. 1, a rough idea of the practical value of Az can 
be obtained. The symptoms which result from too large a Az selection are, 
obviously, failure to obtain a convergent solution or the failure to even 
start the iteration process. 
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At what point does the heating prediction obtained by the SlA Code 
constitute a "solution" to a particular entry problem? This question is at 
the heart of the convergence difficulties. For nonblowing problems the 
author has normally set a convergence level of 5$ with the realization that 
the precision in the radiative and convective flux calculation will be at least 
that good. For a blowing problem the author has worked toward a goal of a 
10$ level of precision. By this is meant the following: Let the radiative 

flux at some cycle be the i^ h value. Then, if the computation were continued 
an endless number of iterations, the asymptotically approached radiative 
heating value would be within 10$ of the i^ value. The degree of precision 
is measured by comparing the radiative flux from the guessed and calculated 
enthalpy. The user is warned, however, that in order for the above measured 
precision to be valid the enthalpy profile must be converged to roughly 10$ 
at most points across the shock layer. The exception is the few points in the 
shear layer region where large spatial variations occur: In the final analysis 

it is the responsibility of the user to interpret the results of the SL^ 

Code in an intelligent way. 
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